source('helpers.R')
library(MaxPac)
library(effectsize)
theme_set(theme_minimal())


suppressWarnings(dir.create('PaperFigures'))

Reading data

We read data from the RetroDuration questionnaire. Translate responses using QTranslate() (google sheet). Read response to the initial_RetrospectiveDuration question. Extract estimated duration since last login as ##:##:## –> Hours:Minutes:Seconds. Compute ClockDuration as the difference between current time and last login.

Compute :

load(file.path(dirBlursday,'Daily_Login_Times.RData'))

# Extracting data:
# 
# RetroDuration <- gimmeRdata(dirBlursday, UniqueName = 'RetroDuration') %>%
#   filter(!Country %in% c('UK','CO', 'US')) %>%
#   QTranslate() %>%
#   filter(Question_Key %in% c('initial_RetrospectiveDuration')) %>%
#   mutate(Date = lubridate::date(Local_Date)) %>%
#   pivot_wider(id_cols = c('PID', 'Session','Run', 'Country', 'Local_Date'), names_from = Question_Key, values_from = Response) %>%
#   # take only the first response if there are several ones
#   group_by(PID,Session,Run,Country) %>%
#   slice(1) %>%
#   add_Demographics() %>%
#   # For EstimatedDuration
#   # Using 3 AM to change day
#   mutate(Date = lubridate::as_date(Local_Date - lubridate::duration(hours = 3))) %>%
#   left_join(select(daily_login_times, PID, Session, Local_Date, Date), by = c('PID','Date'), suffix = c('','_First_Login')) %>%
#   mutate(Date = lubridate::as_date(Local_Date)) %>%
#   mutate(ClockDuration = Local_Date - Local_Date_First_Login) %>%
#   select(Country, PID, Session, Run, Local_Date, ClockDuration, initial_RetrospectiveDuration,
#          Handedness, Sex, Age, Date) %>%
#   mutate(M = str_match(initial_RetrospectiveDuration,
#                        '^.*?(\\d+).*?[:：].*?(\\d+).*?[:：].*?(\\d+).*?$'),
#          h = suppressWarnings(as.numeric(M[,2])),
#          m = suppressWarnings(as.numeric(M[,3])),
#          s = suppressWarnings(as.numeric(M[,4])),
#          EstimatedDuration = lubridate::duration(hour = h, minute = m, second = s),
#          ClockDuration = lubridate::as.duration(ClockDuration),
#          rEstimatedDuration = EstimatedDuration / ClockDuration,
#          EstimationError = EstimatedDuration - ClockDuration,
#          rEstimationError = EstimationError / ClockDuration) %>%
#   select(-M) %>% # ,-h,-m,-s
#   ungroup()
# 
# save(RetroDuration, file = file.path(dirBlursday,'RetroDuration.RData'))

load(file = file.path(dirBlursday,'RetroDuration.RData'))

Outliers detection and clean up

RetroDuration %>%
  mutate(across(c(ClockDuration, EstimatedDuration, rEstimatedDuration, rEstimationError), as.numeric)) %>%
  pivot_longer(cols = c(ClockDuration, EstimatedDuration, rEstimatedDuration, rEstimationError)) %>%
ggplot(aes(x = value, fill = Country)) + 
  geom_histogram(bins = 50) +
  facet_wrap(~name ,scales = 'free') +
  ggtitle('Measure distributions before outliers detection')
## Warning: Removed 1464 rows containing non-finite values (stat_bin).

Distributions of ClockDuration, EstimatedDuration, rEstimationError, and rEstimatedDuration are highly skewed to the left. First capping EstimatedDuration and ClockDuration to a decent duration, then removing 5% most extreme rEstimationError.

Rationale:

We count these

# capped_outliers:
Duration_max <- 18000
Duration_min <- 60
Ratio_rejection_ClockDuration_to_EstimatedDuration <- 5
# trim_outliers:
probs <- c(.025, .975)


# RetroDuration_clean <- RetroDuration %>%
#   filter(!is.na(rEstimatedDuration), rEstimatedDuration != 0)
# nrowsbefore_outliers <- nrow(RetroDuration_clean)
# RetroDuration_clean <- RetroDuration_clean %>%
#   group_by(Country,Session) %>%
#   mutate(outliers = ClockDuration > Duration_max | ClockDuration < Duration_min | EstimatedDuration > Duration_max * Ratio_rejection_ClockDuration_to_EstimatedDuration | EstimatedDuration < Duration_min / Ratio_rejection_ClockDuration_to_EstimatedDuration) %>%
#   filter(!outliers) %>%
#   mutate(outliers = find.outliers(rEstimationError,meth = 'trim', probs = probs)) %>%
#   filter(!outliers) %>%
#   ungroup() %>% select(-outliers)
# nrowsafter_outliers <- nrow(RetroDuration_clean)
# 
# cat('Global proportion of outliers', prop_outliers_global <- 1 - nrowsafter_outliers / nrowsbefore_outliers)
# 
# save(RetroDuration_clean, file = file.path(dirBlursday,'RetroDuration_clean.RData'))
load(file.path(dirBlursday,'RetroDuration_clean.RData'))

RetroDuration_clean %>%
  group_by(Country,Session) %>%
  mutate(across(c(ClockDuration, EstimatedDuration, rEstimatedDuration, rEstimationError), as.numeric)) %>%
  pivot_longer(cols = c(ClockDuration, EstimatedDuration, rEstimatedDuration, rEstimationError)) %>%
ggplot(aes(x = value, fill = Country)) + 
  geom_histogram(bins = 50) +
  facet_wrap(~name ,scales = 'free') +
  ggtitle('Measure distributions after outliers detection')

Subjects Count after outlier rejection

RetroDuration_clean %>% group_by(Country, Session) %>%
  summarize(n = n())

ratio rEstimatedDuration = EstimatedDuration / ClockDuration

Focus on S1 vs. SC comparison

load(file = file.path(dirBlursday, 'RetroDuration_clean.RData'))

tostat <- RetroDuration_clean %>% 
  group_by(Country, PID, Session) %>%
  select(Country,PID, Session, Local_Date,EstimatedDuration, rEstimatedDuration, EstimationError, Age, Sex, Handedness, ClockDuration) %>%
  filter(Session %in% c('S1','SC')) %>%
  summarize(Local_Date = unique(Local_Date),
            EstimationError = mean(EstimationError) / 60,
            ClockDuration = mean(ClockDuration) / 60,
            EstimatedDuration = mean(EstimatedDuration) / 60,
            rEstimatedDuration = mean(rEstimatedDuration),
            rEstimationError = mean(EstimationError/ClockDuration),
            Age = unique(Age)) %>%
  mutate(logrEstimatedDuration = log10(rEstimatedDuration)) %>%
  add_SubjectiveConfinementIndices() %>%
  add_StringencyIndex() %>%
  add_TimeOfDay()

N <- tostat %>%
  group_by(Session) %>%
  summarize(N = n())

Model log(EstimatedDuration) = f(log(ClockDuration))

m <- lm(log10(EstimatedDuration) ~ log10(ClockDuration) * Session, data = tostat)
plot(m)

summary(m)
## 
## Call:
## lm(formula = log10(EstimatedDuration) ~ log10(ClockDuration) * 
##     Session, data = tostat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.84487 -0.10233  0.01158  0.11975  1.10009 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     0.12709    0.02891   4.396 1.17e-05 ***
## log10(ClockDuration)            0.89451    0.01938  46.159  < 2e-16 ***
## SessionSC                       0.16688    0.06903   2.418   0.0157 *  
## log10(ClockDuration):SessionSC -0.12640    0.04594  -2.751   0.0060 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1804 on 1739 degrees of freedom
## Multiple R-squared:  0.587,  Adjusted R-squared:  0.5863 
## F-statistic:   824 on 3 and 1739 DF,  p-value: < 2.2e-16
anova(m)
effectsize::eta_squared(m)
confint(m)
##                                      2.5 %      97.5 %
## (Intercept)                     0.07038706  0.18379810
## log10(ClockDuration)            0.85650572  0.93252236
## SessionSC                       0.03149689  0.30227052
## log10(ClockDuration):SessionSC -0.21650737 -0.03629085
(fig1a <- emmeans(m,~ClockDuration+Session,at = list(ClockDuration = seq(min(tostat$ClockDuration,na.rm = F), max(tostat$ClockDuration,na.rm = F), length.out = 20))) %>%
   summary(type = 'response') %>%
   ggplot(aes(x = ClockDuration, y = response, ymin = lower.CL, ymax = upper.CL, col = Session, group = Session)) +
   geom_abline(slope = 1) +
   geom_point(aes(y= EstimatedDuration, ymin = NA, ymax = NA),alpha = .5, size = 1, data = tostat) +
   geom_line(size = 1) +
   geom_ribbon(col = NA, fill = 'grey',alpha = .2) +
   scale_x_log10() +
   scale_y_log10() +
   ylab('Duration Estimate (min)') +
   xlab('Clock Duration (min)') +
   scale_fill_manual(values = Palette_SessionS1SC) +
   scale_color_manual(values = Palette_SessionS1SC))
## Warning: Ignoring unknown aesthetics: ymin, ymax

summary(m)
## 
## Call:
## lm(formula = log10(EstimatedDuration) ~ log10(ClockDuration) * 
##     Session, data = tostat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.84487 -0.10233  0.01158  0.11975  1.10009 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     0.12709    0.02891   4.396 1.17e-05 ***
## log10(ClockDuration)            0.89451    0.01938  46.159  < 2e-16 ***
## SessionSC                       0.16688    0.06903   2.418   0.0157 *  
## log10(ClockDuration):SessionSC -0.12640    0.04594  -2.751   0.0060 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1804 on 1739 degrees of freedom
## Multiple R-squared:  0.587,  Adjusted R-squared:  0.5863 
## F-statistic:   824 on 3 and 1739 DF,  p-value: < 2.2e-16
summary(emtrends(m,~Session, 'log10(ClockDuration)', infer = T, offset = -1, adjust = 'tukey'))
## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons
## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons
anova(m)
qqnorm(residuals(m))

effectsize::eta_squared(m)
summary(emmeans(m,~ClockDuration+Session,at = list(ClockDuration = c(10,30,60))),type = 'response') %>%
  mutate(ClockDuration = lubridate::as.duration(ClockDuration * 60),
         SE = SE * 60,
         response = lubridate::as.duration(response * 60))
summary(emmeans(m,~ClockDuration+Session,at = list(ClockDuration = seq(15,20, by = .1))),type = 'response') %>%
  mutate(ClockDuration = lubridate::as.duration(ClockDuration * 60),
         response = lubridate::as.duration(response * 60),
         identity = response - ClockDuration) %>%
  group_by(Session) %>%
  filter(abs(identity) == min(abs(identity))) %>%
  select(-response)
summary(emmeans(m,~ClockDuration, at = list(ClockDuration = seq(15,20, by = .1))),type = 'response') %>%
  mutate(ClockDuration = lubridate::as.duration(ClockDuration * 60),
         response = lubridate::as.duration(response * 60),
         identity = response - ClockDuration) %>%
  filter(abs(identity) == min(abs(identity))) %>%
  select(-response)
## NOTE: Results may be misleading due to involvement in interactions
tostat %>% group_by(Session) %>% summarize(mean_ClockDuration = mean(ClockDuration))

Effect of age and Hour of day

tostat %>%
  ggplot(aes(x=Age,y = after_stat(ncount), fill = Session)) +
  geom_histogram(position = 'dodge') +
  scale_fill_manual(values = Palette_SessionS1SC) +
  scale_color_manual(values = Palette_SessionS1SC) +
  ylab('normalized count')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 45 rows containing non-finite values (stat_bin).

All data and sessions

We first insert all covariates in the data, looking at all sessions, now.

load(file = file.path(dirBlursday, 'RetroDuration_clean.RData'))

tostat <- RetroDuration_clean %>% 
  group_by(Country, PID, Session) %>%
  select(Country,PID, Session, Local_Date, EstimatedDuration, rEstimatedDuration, EstimationError, Age, Sex, Handedness, ClockDuration) %>%
  filter(!is.na(rEstimatedDuration), rEstimatedDuration != 0) %>%
  summarize(n = n(),
            Local_Date = unique(Local_Date),
            EstimatedDuration = mean(EstimatedDuration),
            rEstimatedDuration = mean(rEstimatedDuration),
            rEstimationError = mean(EstimationError/ClockDuration),
            EstimationError = mean(EstimationError),
            ClockDuration = mean(ClockDuration),
            Age = unique(Age)) %>%
  mutate(logrEstimatedDuration = log10(rEstimatedDuration)) %>%
  add_SubjectiveConfinementDuration() %>%
  add_SubjectiveConfinementIndices() %>%
  add_Mobility() %>%
  add_StringencyIndex() %>%
  add_TimeOfDay()
## Joining, by = c("Country", "PID")
(N <- tostat %>%
  group_by(Session) %>%
  summarize(N = n()))

Modeling

Subjects were repeatedly tested in S1-S4, and a different pool of subjects were tested in SC. We thus try to see if the error related to subjects in the intercept of EstimatedDuration is worth modeling as a random effect. The standard deviation of the estimated random effects is much smaller than that of the residuals, suggesting that random effects can be ignored in this case.

m <- lmerTest::lmer(log10(rEstimatedDuration) ~ Stringency_Index + Subjective_Confinement + Age + poly(Hour_Of_Day,3) +  (1 | PID), data = tostat )
summary(m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## log10(rEstimatedDuration) ~ Stringency_Index + Subjective_Confinement +  
##     Age + poly(Hour_Of_Day, 3) + (1 | PID)
##    Data: tostat
## 
## REML criterion at convergence: -843.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.5107 -0.5250  0.0382  0.6064  8.3086 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PID      (Intercept) 0.003692 0.06076 
##  Residual             0.035068 0.18726 
## Number of obs: 2147, groups:  PID, 1736
## 
## Fixed effects:
##                          Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)             3.513e-02  2.882e-02  1.853e+03   1.219  0.22304   
## Stringency_Index       -6.230e-04  2.310e-04  1.699e+03  -2.697  0.00706 **
## Subjective_Confinement -7.860e-04  1.309e-03  1.935e+03  -0.600  0.54829   
## Age                     5.869e-05  2.802e-04  1.597e+03   0.209  0.83413   
## poly(Hour_Of_Day, 3)1  -1.730e-01  2.089e-01  2.126e+03  -0.828  0.40760   
## poly(Hour_Of_Day, 3)2  -3.149e-01  2.055e-01  2.139e+03  -1.532  0.12560   
## poly(Hour_Of_Day, 3)3   3.172e-01  2.070e-01  2.139e+03   1.533  0.12553   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Strn_I Sbjc_C Age    p(H_O_D,3)1 p(H_O_D,3)2
## Strngncy_In -0.662                                             
## Sbjctv_Cnfn -0.691  0.088                                      
## Age         -0.361  0.017 -0.027                               
## p(H_O_D,3)1 -0.048  0.014  0.024  0.059                        
## p(H_O_D,3)2 -0.061  0.065  0.028  0.004  0.011                 
## p(H_O_D,3)3 -0.003  0.051 -0.014 -0.044 -0.040       0.011
anova(m)
effectsize::eta_squared(m)

Indeed, running a simple lm on the data with the same fixed effects retrieves the same effects.

m <- lm(log10(rEstimatedDuration) ~ Stringency_Index + Subjective_Confinement + Age + poly(Hour_Of_Day,3), data = tostat)
summary(m)
## 
## Call:
## lm(formula = log10(rEstimatedDuration) ~ Stringency_Index + Subjective_Confinement + 
##     Age + poly(Hour_Of_Day, 3), data = tostat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2366 -0.1115  0.0067  0.1270  1.7101 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)   
## (Intercept)             4.374e-02  2.835e-02   1.543  0.12300   
## Stringency_Index       -6.727e-04  2.259e-04  -2.978  0.00293 **
## Subjective_Confinement -9.998e-04  1.292e-03  -0.774  0.43919   
## Age                     2.528e-05  2.734e-04   0.092  0.92634   
## poly(Hour_Of_Day, 3)1  -1.439e-01  2.083e-01  -0.691  0.48972   
## poly(Hour_Of_Day, 3)2  -2.736e-01  2.055e-01  -1.332  0.18307   
## poly(Hour_Of_Day, 3)3   3.098e-01  2.070e-01   1.497  0.13463   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.197 on 2140 degrees of freedom
##   (267 observations deleted due to missingness)
## Multiple R-squared:  0.006255,   Adjusted R-squared:  0.003469 
## F-statistic: 2.245 on 6 and 2140 DF,  p-value: 0.03654
anova(m)
effectsize::eta_squared(m)
plot(m)

Same reasoning with Mobility_Transit:

m <- lmerTest::lmer(log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit + Subjective_Confinement + Age  +  poly(Hour_Of_Day,3) + (1 | PID), data = tostat )
## Warning: Some predictor variables are on very different scales: consider
## rescaling

## Warning: Some predictor variables are on very different scales: consider
## rescaling
summary(m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit +  
##     Subjective_Confinement + Age + poly(Hour_Of_Day, 3) + (1 |      PID)
##    Data: tostat
## 
## REML criterion at convergence: -836.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.3938 -0.5133  0.0453  0.6010  8.4372 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PID      (Intercept) 0.003564 0.0597  
##  Residual             0.035079 0.1873  
## Number of obs: 2147, groups:  PID, 1736
## 
## Fixed effects:
##                          Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)             5.536e-02  2.972e-02  1.685e+03   1.863 0.062654 .  
## Stringency_Index       -1.408e-03  3.726e-04  1.417e+03  -3.780 0.000163 ***
## Mobility_Transit       -7.562e-04  2.824e-04  2.091e+03  -2.678 0.007464 ** 
## Subjective_Confinement -7.120e-04  1.307e-03  1.932e+03  -0.545 0.585982    
## Age                    -1.624e-05  2.810e-04  1.584e+03  -0.058 0.953914    
## poly(Hour_Of_Day, 3)1  -1.815e-01  2.086e-01  2.125e+03  -0.870 0.384345    
## poly(Hour_Of_Day, 3)2  -3.008e-01  2.052e-01  2.138e+03  -1.466 0.142868    
## poly(Hour_Of_Day, 3)3   3.051e-01  2.067e-01  2.138e+03   1.476 0.140090    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Strn_I Mblt_T Sbjc_C Age    p(H_O_D,3)1 p(H_O_D,3)2
## Strngncy_In -0.594                                                    
## Mblty_Trnst -0.251  0.786                                             
## Sbjctv_Cnfn -0.663  0.037 -0.023                                      
## Age         -0.373  0.088  0.098 -0.029                               
## p(H_O_D,3)1 -0.050  0.022  0.017  0.023  0.060                        
## p(H_O_D,3)2 -0.053  0.022 -0.023  0.029  0.002  0.011                 
## p(H_O_D,3)3 -0.008  0.049  0.021 -0.015 -0.042 -0.039       0.011     
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
anova(m)
effectsize::eta_squared(m)
## Warning: Some predictor variables are on very different scales: consider
## rescaling

Indeed, running a simple lm on the data with the same fixed effects retrieves the same effects.

m <- lm(log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit + Subjective_Confinement + Age + poly(Hour_Of_Day,3), data = tostat)
summary(m)
## 
## Call:
## lm(formula = log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit + 
##     Subjective_Confinement + Age + poly(Hour_Of_Day, 3), data = tostat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.21370 -0.10848  0.00998  0.12468  1.73176 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             6.240e-02  2.909e-02   2.145  0.03208 *  
## Stringency_Index       -1.457e-03  3.616e-04  -4.029  5.8e-05 ***
## Mobility_Transit       -7.790e-04  2.808e-04  -2.774  0.00558 ** 
## Subjective_Confinement -9.164e-04  1.291e-03  -0.710  0.47772    
## Age                    -4.676e-05  2.742e-04  -0.171  0.86461    
## poly(Hour_Of_Day, 3)1  -1.559e-01  2.080e-01  -0.750  0.45354    
## poly(Hour_Of_Day, 3)2  -2.607e-01  2.052e-01  -1.271  0.20402    
## poly(Hour_Of_Day, 3)3   2.979e-01  2.067e-01   1.441  0.14965    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1967 on 2139 degrees of freedom
##   (267 observations deleted due to missingness)
## Multiple R-squared:  0.009818,   Adjusted R-squared:  0.006578 
## F-statistic:  3.03 on 7 and 2139 DF,  p-value: 0.003583
anova(m)
effectsize::eta_squared(m)
confint(m)
##                               2.5 %        97.5 %
## (Intercept)             0.005345989  0.1194542821
## Stringency_Index       -0.002166144 -0.0007477780
## Mobility_Transit       -0.001329668 -0.0002283539
## Subjective_Confinement -0.003447168  0.0016144056
## Age                    -0.000584506  0.0004909846
## poly(Hour_Of_Day, 3)1  -0.563851130  0.2519801769
## poly(Hour_Of_Day, 3)2  -0.663098509  0.1416838641
## poly(Hour_Of_Day, 3)3  -0.107437203  0.7032299584
plot(m)

### Illustrating

Stringency Index

m <- lm(log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit + Subjective_Confinement + Age + poly(Hour_Of_Day,3), data = tostat)

em <- summary(emmeans(m, ~ Stringency_Index, at = list(Stringency_Index = seq(20,100,by = 20))), type = 'response')

(fig1b <- ggplot(aes(x = Stringency_Index,y = rEstimatedDuration, col = Session), data = tostat) +
  geom_jitter(alpha = .3, size = 1) +
  geom_line(aes(y=response), col = 'black', data = em) +
  geom_ribbon(aes(y=response, ymin = lower.CL, ymax = upper.CL), col = NA, alpha = .2, data = em) +
  scale_y_log10( limits = c(.3,3)) +
  scale_fill_manual(values = Palette_Session) +
  scale_color_manual(values = Palette_Session) +
  ylab('Relative Duration Estimate') +
  xlab('Stringency Index'))
## Warning: Removed 36 rows containing missing values (geom_point).

#### Mobility Transit

em <- summary(emmeans(m, ~ Mobility_Transit, at = list(Mobility_Transit = seq(-100,25,by = 10))), type = 'response')

(fig1c <- ggplot(aes(x = Mobility_Transit,y = rEstimatedDuration, col = Session), data = tostat) +
  geom_jitter(alpha = .3, size = 1) +
  geom_line(aes(y=response), col = 'black', data = em) +
  geom_ribbon(aes(y=response, ymin = lower.CL, ymax = upper.CL), col = NA, alpha = .2, data = em) +
  scale_y_log10( limits = c(.3,3)) +
  scale_fill_manual(values = Palette_Session) +
  scale_color_manual(values = Palette_Session) +
  ylab('Relative Duration Estimate') +
  xlab('Mobility Index'))
## Warning: Removed 36 rows containing missing values (geom_point).

Other mobility measures have little effect on rEstimatedDuration

m <- lm(log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit + Mobility_Retail + Mobility_Parks + Mobility_WorkPlaces + Mobility_Residential + Subjective_Confinement + Age + poly(Hour_Of_Day,3), data = tostat)
summary(m)
## 
## Call:
## lm(formula = log10(rEstimatedDuration) ~ Stringency_Index + Mobility_Transit + 
##     Mobility_Retail + Mobility_Parks + Mobility_WorkPlaces + 
##     Mobility_Residential + Subjective_Confinement + Age + poly(Hour_Of_Day, 
##     3), data = tostat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.20132 -0.10766  0.00957  0.12317  1.72235 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)  
## (Intercept)             5.222e-02  3.418e-02   1.528   0.1267  
## Stringency_Index       -9.978e-04  4.844e-04  -2.060   0.0395 *
## Mobility_Transit       -1.331e-03  7.323e-04  -1.817   0.0694 .
## Mobility_Retail         1.069e-03  7.004e-04   1.527   0.1270  
## Mobility_Parks         -3.610e-05  2.142e-04  -0.169   0.8661  
## Mobility_WorkPlaces     1.122e-03  8.110e-04   1.383   0.1667  
## Mobility_Residential    3.119e-03  1.845e-03   1.690   0.0911 .
## Subjective_Confinement -1.217e-03  1.306e-03  -0.931   0.3518  
## Age                    -8.432e-05  2.776e-04  -0.304   0.7613  
## poly(Hour_Of_Day, 3)1  -1.612e-01  2.084e-01  -0.774   0.4392  
## poly(Hour_Of_Day, 3)2  -2.783e-01  2.055e-01  -1.355   0.1757  
## poly(Hour_Of_Day, 3)3   3.045e-01  2.076e-01   1.467   0.1425  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1967 on 2135 degrees of freedom
##   (267 observations deleted due to missingness)
## Multiple R-squared:  0.01226,    Adjusted R-squared:  0.007167 
## F-statistic: 2.408 on 11 and 2135 DF,  p-value: 0.005686
anova(m)
effectsize::eta_squared(m)
plot(m)

Time Of Day

Figure1 final

library(patchwork)
## 
## Attaching package: 'patchwork'
## The following object is masked from 'package:cowplot':
## 
##     align_plots
layout <- "
A
B
B"

fig1a / (fig1b / fig1c + plot_layout(guides = 'collect')) + plot_layout(design = layout) + plot_annotation(tag_levels = 'a') & theme(plot.tag = element_text(face = 'bold'))
## Warning: Removed 36 rows containing missing values (geom_point).
## Removed 36 rows containing missing values (geom_point).

imagefile = 'PaperFigures/Fig1.pdf'
ggsave(
  basename(imagefile),
  plot = last_plot(),
  device = 'pdf',
  width = 89,
  height = 150,
  path = dirname(imagefile),
  units = 'mm',
  bg = "white"
)
## Warning: Removed 36 rows containing missing values (geom_point).
## Removed 36 rows containing missing values (geom_point).